Spatially Continuous Data IV

Spatially Continuous Data IV

NOTE: You can download the source files for this book from here. The source files are in the format of R Notebooks. Notebooks are pretty neat, because the allow you execute code within the notebook, so that you can work interactively with the notes.

In the previous practice you were introduced to the concept of variographic analysis for fields/spatially continuous data.

If you wish to work interactively with this chapter you will need the following:

• An R markdown notebook version of this document (the source file).

• A package called geog4ga3.

Learning objectives

In this practice, you will learn:

1. Using residual spatial pattern to estimate prediction errors.
2. Kriging: a method for optimal predictions.

Suggested reading

• Bailey TC and Gatrell AC [-@Bailey1995] Interactive Spatial Data Analysis, Chapters 5 and 6. Longman: Essex.
• Bivand RS, Pebesma E, and Gomez-Rubio V [-@Bivand2008] Applied Spatial Data Analysis with R, Chapter 8. Springer: New York.
• Brunsdon C and Comber L [-@Brunsdon2015R] An Introduction to R for Spatial Analysis and Mapping, Chapter 6, Sections 6.7 and 6.8. Sage: Los Angeles.
• Isaaks EH and Srivastava RM [-@Isaaks1989applied] An Introduction to Applied Geostatistics, CHAPTERS. Oxford University Press: Oxford.
• O’Sullivan D and Unwin D [-@Osullivan2010] Geographic Information Analysis, 2nd Edition, Chapters 9 and 10. John Wiley & Sons: New Jersey.

Preliminaries

As usual, it is good practice to clear the working space to make sure that you do not have extraneous items there when you begin your work. The command in R to clear the workspace is rm (for “remove”), followed by a list of items to be removed. To clear the workspace from all objects, do the following:

rm(list = ls())

Note that ls() lists all objects currently on the worspace.

Load the libraries you will use in this activity:

library(tidyverse)
library(spdep)
library(plotly)

Attaching package: 㤼㸱plotly㤼㸲

The following object is masked from 㤼㸱package:ggplot2㤼㸲:

    last_plot

The following object is masked from 㤼㸱package:stats㤼㸲:

    filter

The following object is masked from 㤼㸱package:graphics㤼㸲:

    layout

library(gstat)
library(geog4ga3)

Begin by loading the data file:

data("Walker_Lake")

You can verify the contents of the dataframe:

summary(Walker_Lake)

      ID                  X               Y               V                U           T      
 Length:470         Min.   :  8.0   Min.   :  8.0   Min.   :   0.0   Min.   :   0.00   1: 45  
 Class :character   1st Qu.: 51.0   1st Qu.: 80.0   1st Qu.: 182.0   1st Qu.:  83.95   2:425  
 Mode  :character   Median : 89.0   Median :139.5   Median : 425.2   Median : 335.00          
                    Mean   :111.1   Mean   :141.3   Mean   : 435.4   Mean   : 613.27          
                    3rd Qu.:170.0   3rd Qu.:208.0   3rd Qu.: 644.4   3rd Qu.: 883.20          
                    Max.   :251.0   Max.   :291.0   Max.   :1528.1   Max.   :5190.10          
                                                                     NA's   :195              

Using residual spatial pattern to estimate prediction errors

Previously, you have seen how to interpolate a field using trend surface analysis. You have also seen how that may at least on occasions residuals that are not spatially independent.

The implication of non-random residuals is that there uncaptured pattern remains. Which means that information can still be extracted from it. Lets revisit the case of Walker Lake to explore one way of doing this.

As before, we first calculate the polynomial terms of the coordinates to fit a trend surface to the data:

Walker_Lake <- mutate(Walker_Lake,
                        X3 = X^3, X2Y = X^2 * Y, X2 = X^2, 
                        XY = X * Y,
                        Y2 = Y^2, XY2 = X * Y^2, Y3 = Y^3)

Given the polynomial expansion, we can proceed to estimate the following cubic trend surface model, which in the previous practice provided the best fit to the data:

WL.trend3 <- lm(formula = V ~ X3 + X2Y + X2 + X + XY + Y + Y2 + XY2 + Y3, 
                data = Walker_Lake)
summary(WL.trend3)

Call:
lm(formula = V ~ X3 + X2Y + X2 + X + XY + Y + Y2 + XY2 + Y3, 
    data = Walker_Lake)

Residuals:
    Min      1Q  Median      3Q     Max 
-564.19 -197.41    7.91  194.25  929.72 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -8.620e+00  1.227e+02  -0.070 0.944035    
X3           1.533e-04  4.806e-05   3.190 0.001522 ** 
X2Y          6.139e-05  3.909e-05   1.570 0.117000    
X2          -6.651e-02  1.838e-02  -3.618 0.000330 ***
X            9.172e+00  2.386e+00   3.844 0.000138 ***
XY          -4.420e-02  1.430e-02  -3.092 0.002110 ** 
Y            4.794e+00  2.040e+00   2.350 0.019220 *  
Y2          -1.806e-03  1.327e-02  -0.136 0.891822    
XY2          7.679e-05  2.956e-05   2.598 0.009669 ** 
Y3          -4.170e-05  2.819e-05  -1.479 0.139759    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 276.7 on 460 degrees of freedom
Multiple R-squared:  0.1719,    Adjusted R-squared:  0.1557 
F-statistic: 10.61 on 9 and 460 DF,  p-value: 5.381e-15

Lets visualize the residuals. As you can see, the residuals do not appear to be random

plot_ly(x = ~Walker_Lake$X, y = ~Walker_Lake$Y, z = ~WL.trend3$residuals, color = ~WL.trend3$residuals < 0, colors = c("blue", "red"), type = "scatter3d")

No scatter3d mode specifed:
  Setting the mode to markers
  Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
No scatter3d mode specifed:
  Setting the mode to markers
  Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

Now create an interpolation grid.

X.p <- seq(from = 0.1, to = 255.1, by = 2.5)
Y.p <- seq(from = 0.1, to = 295.1, by = 2.5)
df.p <- expand.grid(X = X.p, Y = Y.p)

to which we can add the polynomial terms:

df.p <- mutate(df.p, X3 = X^3, X2Y = X^2 * Y, X2 = X^2, 
               XY = X * Y, 
               Y2 = Y^2, XY2 = X * Y^2, Y3 = Y^3)

The interpolated cubic surface is obtained as:

WL.preds3 <- predict(WL.trend3, newdata = df.p, se.fit = TRUE, interval = "prediction", level = 0.95)

The surface is converted into a matrix for 3D plotting:

z.p3 <- matrix(data = WL.preds3$fit[,1], nrow = length(Y.p), ncol = length(X.p), byrow = TRUE)

And plot:

WL.plot3 <- plot_ly(x = ~X.p, y = ~Y.p, z = ~z.p3, 
        type = "surface", colors = "YlOrRd") %>% 
  layout(scene = list(
    aspectmode = "manual", aspectratio = list(x = 1, y = 1, z = 1)))
WL.plot3

The trend surface provides a smooth estimate of the field. However, it is not sufficient to capture all systematic variation, and fails to produce random residuals.

A possible way of enhancing this approach to interpolation is to exploit the information that remains in the residuals, for instance by the use of \(k\)-point means.

Lets see how this is done.

To interpolate the residuals, we first need the set of target points (the points for the interpolation), as well as the source (the observations)m as follows:

target_xy = expand.grid(x = X.p, y = Y.p)
source_xy = cbind(x = Walker_Lake$X, y = Walker_Lake$Y)

It is possible now to use the kpointmean function to interpolate the residuals, for instance using \(k=5\) neighbors:

kpoint.5 <- kpointmean(source_xy = source_xy, z = WL.trend3$residuals, target_xy = target_xy, k = 5)

Given the estimated residuals, we can add them to the cubic trend surface, as follows:

z.p3 <- matrix(data = WL.preds3$fit[,1] + kpoint.5$z,
               nrow = length(Y.p), ncol = length(X.p), byrow = TRUE)

to produce an interpolated field:

WL.plot3 <- plot_ly(x = ~X.p, y = ~Y.p, z = ~z.p3, 
        type = "surface", colors = "YlOrRd") %>% 
  layout(scene = list(
    aspectmode = "manual", aspectratio = list(x = 1, y = 1, z = 1)))
WL.plot3

Of all the approaches you have seen so far, this is the first that provides a genuine estimate of the following: \[ \hat{z}_p + \hat{\epsilon}_p \]

With trend surface analysis providing a smooth estimator of the underlying field: \[ \hat{z}_p = f(x_p, y_p) \]

And \(k\)-point means providing an estimator of: \[ \hat{\epsilon}_p \]

A question is how to decide the number of neighbors to use in the calculation of the k-point means. As previously discussed, \(k = 1\) becomes identical to Voronoi polygons, and \(k = n\) becomes the global mean.

A second question concerns the way the average is calculated. As variographic analysis demonstrates, it is possible to estimate the way in which spatial dependence weakens with distance. Why should more distant points be weighted equally? The answer is, there is no reason why they should, and in fact, variographic analysis elegantly solves this, as well the question of how many points to use: all of them, with varying weights.

Next, we will introduce kriging, a method for optimal prediction that is based on the use of variographic analyisis.

Kriging: a method for optimal prediction.

Lets begin by positing a situation as follows: \[ \hat{z}_p + \hat{\epsilon}_p = \hat{f}(x_p, y_p) + \hat{\epsilon}_p \]

where \(\hat{f}(x_p, y_p)\) is a smooth estimator of an underlying field.

We aim to predict \(\hat{\epsilon}_p\) based on the observed residuals. We use an expression similar to the one used for IDW and k-point means (we will use \(\lambda\) for the weights to avoid confusing the the weights in variographic analysis): \[ \hat{\epsilon}_p = \sum_{i=1}^n {\lambda_{pi}\epsilon_i} \]

That is, \(\hat{\epsilon}_p\) is a linear combination of the prediction residuals from the trend: \[ \epsilon_i = z_i - \hat{f}(x_i, y_i) \]

Since the residuals are random, the prediction residuals are random too, and it is possible to define the following expected mean squared error, or prediction variance: \[ \sigma_{\epsilon}^2 = E[(\hat{\epsilon}_p - \epsilon_i)^2] \]

The prediction variance measures how close, on average, the prediction error is to the residuals.

The prediction variance can be decomposed as follows: \[ \sigma_{\epsilon}^2 = E[\hat{\epsilon}_p] + E[\epsilon_i] - 2E[\hat{\epsilon_i\epsilon}_p] \]

It turns out (we will not show de detailed derivation, but it can be consulted here), that the expresion for the prediction variance depends on the weights: \[ \sigma_{\epsilon}^2 = \sum_{i=1}^n \sum_{j=1}^n{\lambda_{ip}\lambda_{jp}C_{ij}} + \sigma^2 + 2\sum_{i=1}^{n}{\lambda_{ip}C_{ip}} \] where \(C_{ij}\) is the autocovariance between observations at \(i\) and \(j\), and \(C_{ip}\) is the autocovariance between the observation at \(i\) and prediction location \(p\).

Fortunately for us, the semivariogram and the autocovariance is straightforward: \[ C_{z}(h) =\sigma^2 - \hat{\gamma}_{z}(h) \]

This means that, given the distance \(h\) between \(i\) and \(j\), and \(i\) and \(p\), we can use a semivariogram to obtain the autocovariances needed to calculate the prediction variance. We are still missing, however, the weights \(\lambda\), which are not known a priori.

These weights can be obtained if we use the following rules:

The expectation of the prediction errors is zero (unbiassedness) Find the weights \(lambda\) that minimize the prediction variance (optimal estimator).

This makes sense, since we would like our predictions to be unbiased and as precise as possible, that is, to have the smallest variance.

Again, solving the minimization problem is beyond the scope of this practice, but it suffices to say that the results are as follows: \[ \mathbf{\lambda}_p = \mathbf{C}^{-1}\mathit{c}_{p} \] where \(\mathbf{C}\) is the covariance matrix, and \(\mathit{c}_{p}\) is the covariance vector for location \(p\).

Kriging is known to have the properties of Best (in the sense that it minimizes the variance) Linear (because of predictions are a linear combination of weights) Unbiased (since the estimators of the prediction errors are zero) Estimator, or BLUP.

Kriging is implemented in the package gstat as follows.

Lets first conduct variographic analysis of the residuals. The function variogram uses as an argument a SpatialPointsDataFrame that we can create as follows, based on the tidy table:

Walker_Lake.sp <- Walker_Lake
coordinates(Walker_Lake.sp) <- ~X+Y

The variogram of the residuals can be obtained by specifying a trend surface in the formula:

variogram_v <- variogram(V ~ X3 + X2Y + X2 + X + XY + Y + Y2 + XY2 + Y3, 
                         data = Walker_Lake.sp)
ggplot(data = variogram_v, aes(x = dist, y = gamma)) +
  geom_point() + 
  geom_text(aes(label = np), nudge_y = -1500) +
  xlab("Distance") + ylab("Semivariance")

You can verify that the semivariogram above corresponds to the residuals by repeating the analysis directly on the residuals. First join the residuals to the SpatialPointsDataFrame:

Walker_Lake.sp$e <- WL.trend3$residuals

And then calculate the semivariogram and plot:

variogram_e <- variogram(e ~ 1, 
                         data = Walker_Lake.sp)
ggplot(data = variogram_e, aes(x = dist, y = gamma)) +
  geom_point() + 
  geom_text(aes(label = np), nudge_y = -1500) +
  xlab("Distance") + ylab("Semivariance")

The empirical semivariogram is used to estimate a semivariogram function:

variogram_v.t <- fit.variogram(variogram_v, model = vgm("Exp", "Sph", "Gau"))
variogram_v.t

The variogram function plots as follows:

gamma.t <- variogramLine(variogram_v.t, maxdist = 130)
ggplot(data = variogram_v, aes(x = dist, y = gamma)) +
  geom_point(size = 3) + 
  geom_line(data = gamma.t, aes(x = dist, y = gamma)) +
  xlab("Distance") + ylab("Semivariance")

df.sp <- df.p
coordinates(df.sp) <- ~X+Y

V.kriged <- krige(V ~ X3 + X2Y + X2 + X + XY + Y + Y2 + XY2 + Y3,
              Walker_Lake.sp, df.sp, variogram_v.t)

[using universal kriging]

Extract the predictions and prediction variance from the object V.kriged:

V.km <- matrix(data = V.kriged$var1.pred,
               nrow = 119, ncol = 103, byrow = TRUE)
V.sm <- matrix(data = V.kriged$var1.var,
               nrow = 119, ncol = 103, byrow = TRUE)

You can now plot the interpolated field:

V.km.plot <- plot_ly(x = ~X.p, y = ~Y.p, z = ~V.km, 
        type = "surface", colors = "YlOrRd") %>% 
  layout(scene = list(
    aspectmode = "manual", aspectratio = list(x = 1, y = 1, z = 1)))
V.km.plot

Also, you can plot the kriging standard errors (the square root of the prediction variance). This gives an estimate of the uncertainty in the predictions:

V.sm.plot <- plot_ly(x = ~X.p, y = ~Y.p, z = ~sqrt(V.sm), 
        type = "surface", colors = "YlOrRd") %>% 
  layout(scene = list(
    aspectmode = "manual", aspectratio = list(x = 1, y = 1, z = 1)))
V.sm.plot

Where are predictions more/less precise?